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Abstract 

We investigate the critical behaviour of a probabilistic mixture of cellular au- 
tomata (CA) rules 182 and 200 (in Wolfram's enumeration scheme) by mean- 
field analysis and Monte Carlo simulations. We found that as we switch off 
one CA and switch on the other by the variation of the single paramenter of 
the model the probabilistic CA (PCA) goes through an extinction-survival-type 
phase transition, and the numerical data indicate that it belongs to the directed 
percolation universality class of critical behaviour. The PCA displays a char- 
acteristic stationary density profile and a slow, diffusive dynamics close to the 
pure CA 200 point that we discuss briefly. Remarks on an interesting related 
stochastic lattice gas are addressed in the conclusions. 
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1 Introduction 



Cellular automata (CA, also for cellular automaton) have been a major theme of 
investigation since their introduction in the late 40's and early 50's of the last 
century as model systems for parallel processing machines and algorithms, mem- 
ory storage devices, and self -reproducing simple organisms [1]. Meanwhile, CA 
have become deeply intertwined with the subject matter of modern equilibrium and 
nonequilibrium statistical mechanics, such as with the theory of phase transitions, 
irreversibility, ergodicity, chaos, percolation, and physical kinetics, to name a few 
[2,3,4,5,6,7,8,9, 10]. 

Probabilistic versions of CA (PCA) with absorbing states, i.e., with states that 
once reached trap the dynamics definitely, are a whole chapter in CA theory. The 
single most representative PCA in this class is the Domany-Kinzel (DK) PCA, 
because its two-dimensional parameter space encompasses both the directed bond 
and the directed site percolation processes [4] . More general PCA with more than 
one absorbing state have also been investigated, with larger neighbourhoods and 
more than two parameters, at the expense of some symmetries present in the DK 
PCA. In particular, a class of two and three parameter PCA with two absorbing 
states was investigated in [11], where many general results on the phase diagrams 
of the models were obtained. 

In usual PCA modeling, the freedom to set the transition probabilities directly 
in the rule table allows one to embody competing interactions in the PCA from 
the outset, such that a resulting complex dynamic behaviour becomes an expected 
treat. However, elementary one-dimensional CA are more charming in their sim- 
plicity, since most of them can be implemented with a few deterministic binary 
operations that also make them of greater technological relevance. 1 Moreover, it 
has been found that by composing simple CA in space and time, probabilistically 
or not, we may obtain complex behaviour out of simple components. A striking ex- 
ample of this possibility is given by the composite deterministic CA of H. Fuks that 
solves the density classification problem, an impossible problem for locally inter- 
acting single CA [12, 13]. These observations triggered our interest in probabilistic 
mixtures of simple CA. One PCA that we considered preliminarly combines CA 
rules 150 and 200 in Wolfram's enumeration scheme [3] and is briefly mentioned 
in Section 4. Another possible combination is given by CA rules 23 and 254, that 
when combined as a PCA in the guise of a model for an organism (l's) consum- 

! An "atlas" of elementary CA can be consulted at http : / /atlas . wolfram, com/ 01/ 01/. 
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ing finite renewable resources (O's) and dying from overcrowding may display an 
active-inactive type phase transition. As far as we are aware, this kind of composite 
PCA has not yet been fully explored in the literature. 

In this paper we investigate in detail a probabilistic mixture of CA 182 and 
CA 200 that is left-right symmetric, has a single scalar order paramenter, no con- 
served quantities, and two absorbing states, both reachable dynamically as long 
as the mixed PCA contains finite portions of both CA 182 and CA 200 dynam- 
ics. After some exploratory work, we found that this PCA displays some inter- 
esting features, such as its behaviour in a certain small parameter region besides 
an extinction-survival-type phase transition as we probabilistically switch off one 
CA and switch on the other. We characterize this phase transition by mean-field 
analysis and Monte Carlo simulations and determine that it belongs to the directed 
percolation universality class of critical behaviour. 

The article is organized as follows. In Section 2, we define the PCA and in 
Section 3 we analyse it at the mean-field level of approximation. Section 4 presents 
the results of direct Monte Carlo simulations of the PCA for its critical point and 
critical exponents, and in Section 5 we make a few final remarks, summarize our 
results, and identify some perspectives for further investigation. 

2 Model description 

Let rn(t) G {0, 1} denote the state of the site I G A C Z at instant t G N, with A 
a finite lattice of L sites with periodic boundary conditions I + L = I. The state 
of the system at instant t is given by rj(t) = (r/i (t) , r\2 (t) , . . . ,rjL(t)) G {0, 1} A . 
The model we are interested in is the probabilistic cellular automaton defined by 
the rules in Table 1. We dub this system PCA pl82-g200, with p + q = 1. 

When p = 1, the mixed PCA reduces to the deterministic CA 182, which has 
a stationary density of l's given by pg| ct = 3/4 [3]. CA 200 (p = in Table 1), 
otherwise, is quite an uninteresting CA, since most initial configurations die out 
quickly or do not evolve at all under its rules. The role of CA 200 in our mixture 
of CAs is to provide a route to the absorbing state devoid of l's, a state that cannot 
be reached by CA 182 except from the initial empty configuration itself. So, as p 
varies between p = and p = 1 we expect a competition to set up between the two 
rules to dominate the dynamics, with an extinction-survival-type phase transition 
somewhere in between. We indeed found such a phase transition, as long as the 
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Table 1: Rule table for PCA pl82-g200, p + q = 1. The first row gives the initial 
neighborhood, the other two rows give the final state reached by the central bit of the initial 
neighborhood with the probability given at the leftmost column. Clearly, the configurations 
00- ■ -0 and 11- • 1 are absorbing configurations of the PCA. 
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initial condition is not one of the absorbing configurations of the model. 

PCA y>182-g200 can be related with the DK PCA only at p = 1/2, a point that 
is off the mixed site-bond percolation parameter subspace of the DK PCA, since 
this subspace requires that a certain parameter x of the DK PCA be zero, which is 
impossible in our setting, because the match between the two PCAs requires that 
p = 1/2 and x = 1 — p simultaneously. The point p = 1/2 in PCA pl82-g200, 
however, corresponds to a completely uncorrected dynamics, since bits will flip 
with probability 1/2 irrespective of their neighborhood, with the exception of the 
bits in the bulk of blocks 11- • -1, cf. Table 1. These blocks, however, are not stable 
at p ^ 0, since they will eventually be eroded from the boundaries. We thus expect 
that atp = 1/2 the stationary density of active sites fluctuates around p = 1/2. 



3 Mean-field analysis 

We begin by analysing PCA pl82-g200 in the mean-field approximation to obtain 
some first information on its critical behaviour. For a nice brief exposition and 
application of the technique see, e.g., [14]. 

The dynamics of the probability distribution Pt{vi) of the states r/ of the PCA 
is ruled by the equation 

p t+ i(v') = Y< w ( r i'\v)Pt(v), (i) 

where the summation runs over all T) e {0, 1} A and W(rj'\r}) ^ is the condi- 
tional probability for a transition r] — > rj' to occur in one time step. Since in a CA 
or PCA all sites are updated simultaneously and independently, we can write 

L 

W(r t '\v) = l[W e (r, , i \r t ), with ^WeM = 1. (2) 
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For PCA pl82-(7200 we have Wi(rf e \rj) = W(r)'^qt-\, %, independent of 

I. The time evolution of the marginal probability distribution Pt(rj£, i]e+i, 
Ve+n-i) of observing n consecutive sites in state (rj£, . . . , 7#+ n _i) is, from 
Eqs. (1) and (2), given by 

p t+i(Ve,V£+n ■ ■ ■ iVt+n-l) = 

= w (Ve\V£-i,m,Ve+i)W(r)' £+1 \r) £ ,rj£ +1 ,rii + 2) ■ ■ ■ 

1U-1 ru ri e+n 

■ ■ ■ W(r]' e+n _ 1 \r] e+n _2,f]e + n-i,'ni+n) Pt(r}e-i,Vl> ■ ■ -iVl+n)- ( 3 ) 

We see from equation (3) that to determine the probability of observing n con- 
secutive sites in a given state at instant i + lwe need to know the probabilities of 
observing the state of n + 2 sites at instant t. To proceed with the calculations in 
an approximate fashion, we truncate this hierarchy at some point to split the cor- 
relations and get a closed set of equations. The simplest approximation (n = 1) is 
obtained by taking 

P t (ru-x,r]i,Tii+i) « p t{ty-i) Pt{m) Pt(vm)- ( 4 ) 

Higher order approximations (n ^ 2) are obtained (assuming spatial homogeneity) 
by the generalized splitting scheme 

D , \ p t(vi-u ■ ■ ■ ,m+n-2) Ptive, ■ ■ ■ ,m+n-i) Ptive+i, ■ ■ ■ ,ve+n) 

^tW-li • • • >%+nJ 57 7757 v • 

Pt(j)£, ■ ■ -,Ve+n-2) Pt(V£+l: ■ ■ -iW+n-l) 

(5) 

From the rules in Table 1 and the above equations, the single-site (n = 1) approx- 
imation for p^ = Pt(rji = 1) (where the superscript refers to the order of the 
approximation) reads 

P W = (P 4 (1) ) 3 + (2 - P){p?)\l ~ + 3rf } (l - p^f. (6) 
In the stationary state, pf^ = p]/ = pW, and the above equation becomes 

p® = (p«) 3 + (2 - p)(p«) 2 (l ~ + 3pp«(l " P W ) 2 i W 

with solutions pM = 0, pW = 1 and pW = (Sp - l)/(4p - 1). The first 
two solutions correspond to the two absorbing states of the dynamics, whereas 
the last, nontrivial solution corresponds to the active state and is valid as long 
as p ^ 1/3. The single-site mean-field approximation for PCA pl82-g200 thus 
predicts an extinction-survival-type phase transition at p^ = 1/3. Notice that at 
p= 1, pW — 2/3, not far from the exact stationary value p^gf * = 3/4 for CA 182. 
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Figure 1 : Density (p) of active sites in the mean-field approximations of orders n — 
2, 3, and 4, together with Monte Carlo simulation data for a lattice of L = 8000 sites 
initialized randomly with density 1/2. Each symbol in the Monte Carlo curve is an average 
over 10 6 samples. Errors in the data are negligible below the critical point p c ~ 0.48, of 
the order of 5% just at the critical point, and less than 0.2% on the rest of the curve. 



We have also considered higher-order approximations with n = 2, 3, and 4. In 
these cases, however, the equations are too cumbersome to be written down here. 
The n = 2 approximation was solved analytically and gives the same results as the 
n = 1 case — somehow, the two-sites marginal probability P(r]£, rj£ + i) factors into 
P(r]i)P(r]i + i). The cases n = 3 and 4 were solved numerically. The densities 
pip) f ac tive sites for these approximation are shown in Figure 1. As expected, 
at p = 1/2 all approximations (as well as the Monte Carlo simulation results, 
cf. Section 4) give /o(l/2) = 1/2. In Table 2 we show the values of and at 
p = 1 for these approximations together with the corresponding values obtained by 
Monte Carlo simulations. As one can see from Table 2, the n-th order mean-field 
approximation to the exact PCA converges slowly with n. 



4 Direct Monte Carlo simulation 
4.1 The density profile 

Our Monte Carlo simulations of PCA pl82-g200 ran as follows. For each p, the 
PCA is initialized according to a Binomial (L, s) distribution (mostly with s = ^), 
i.e., each site is initialized as rje = 1 with probability s, and relaxed through L 2 /10 
Monte Carlo steps, with one Monte Carlo step equal to a synchronous update of 
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Table 2: Critical parameter p c and value of the density of active sites p^ n \p) at p = 1 
obtained by mean-field approximations or order n and Monte Carlo simulation (cf. Sec- 
tion 4). The numbers between parentheses (in this table and elsewhere in this article) 
indicate the uncertainty in the last digit(s) of the data. For p = 1, PCA pl82-g200 reduces 
to CA 182, for which the exact stationary density Pi|f ct = 3/4. 
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the states of all L sites of the lattice. We then sample pl = L 1 Y^e Vi every other 
Monte Carlo step. 

Our data for Pl(p) in a lattice of L = 8000 sites appear in Figure 2. We can 
identify two regions in this figure: an "active region I" that goes from p = up 
to p ~ 0.03 and an "active region II" that goes from p ~ 0.48 up to p = 1. Re- 
gion I is easily understood as an artifact of our initialization of the PCA with a 
Binomial (L, |) distribution. At p = 0, PCA pl82-g200 becomes CA 200, that 
washes out all local configurations but those of neighbouring l's (cf. Table 1). 
The only neighborhood that evolves under CA 200 dynamics is 010, that becomes 
000 with probability q = 1 — p = 1. For a Binomial (L, s) initial distribu- 
tion, the expected initial density of l's is E[% = 1] = s and of triplets 010 is 
E[r]er]£ + ir]£ + 2 = 010] = (1 — s)s(l — s). After one single time step, all 010 go into 
000 and the stationary density of l's becomes E[r^ = 1] = s — (1 — s)s(l — s). For 
s = 1/2, p(p = 0) = 3/8, in accordance with the value in Figure 2. As p increases 
from zero, the noise brough up by the CA 182 rules disturbs the dynamics and the 
PCA eventually converges to the absorbing state devoid of active sites. 

For very small p > 0, after the triplets 010 quickly become 000 with probability 
q = 1 — p, only the boundaries of the remaining clusters 11- - -1 move, perform- 
ing random walks 100 =± 110 at the right edge and 001 =± 011 at the left edge 
with probability p. Under this dynamics, the number of active sites just fluctuates 
about a certain value. When these wandering boundaries meet, however, either two 
clusters coalesce with a low probability p through the 101 — > 111 channel or one 
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Figure 2: Density (p) of active sites for a lattice of L — 8000 sites initialized randomly 
with density 1/2. Each value of pl(p) is an average over 10 6 samples after relaxation 
through ~ L 2 MCS. Errors in the Monte Carlo data are similar to the ones reported in 
Figure 1. 



cluster vanishes with a high probability 1 — p through the 010 — > 000 channel. 
This behaviour is reminiscent of the long time behaviour of the one-species lattice 
gas where k particles coalesce into £ particles, kA — > IA with k > I, which can 
be mapped into the problem of the reunion of k random walkers bounded to move 
in a limited region [15]. The average density of active sites for very small p is 
then expected to decay very slowly towards zero as p(t) ~ 1/ yjpi or, equivalently, 
the extinction time T ey ±{L,p) ~ L 2 /p, with prefactors depending on the sizes and 
spatial distribution of clusters 11- • -1 on the initial configuration. Since we sam- 
ple only a very tiny fraction of all possible initial configurations (for L = 400 
and p(0) = 0.3, the number of possible configurations under periodic boundary 
conditions is ~ 1.4 x 10 102 ), these prefactors vary widely from one simulation 
to the other, even if we average over several thousands of realizations. This be- 
came evident as the data for T exi (L,p) turned out to be overdispersed, i.e., with 
sample variance exceeding the mean. Numerical evidence for the scaling of r ex t 
is given in Figure 3. We believe that the deviation from the expected exponents is 
due both to finite-size effects and to the prefactor issue. While it is a little both- 
ering to obtain such non-expected figures, at least we got the correct functional 
form T ex t(L,p) ~ L a /p b with a > 1 and 1 » 1, corroborating our qualitative 
analysis. Notice that a ~ zdp — 1.581 (cf. Section 4.2), but since this region of 
small p is very far from the critical region there is no reason to confuse the two 
values — the figures are similar by accident and by the numerical limitations and 
issues mentioned before. 
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Figure 3: Extinction time T ex t(L,p) for small p. Each point shown is an average over 
1000 realizations of the PC A with a random initial distribution of active sites with density 
p(0) = 0.3. The curves fitted are guides only — the prefactors vary widely since only a tiny 
fraction of all possible initial conditions are sampled. 



The unusual bump in pl(p) visible in Figure 2 near p ~ 0.57 does not have 
a simple explanation. We thought at first that it could be signaling a first order 
phase transition, but further simulation data convinced us that the phase transition 
is continuous. Notice that the non-monotonic behaviour of p{p) for p > 0.57 
is also observed in the mean-field approximations of orders n = 3 and 4. We 
found, however, that the closely related PCA pl50-<7200 may provide a clue to 
this bump. PCA pl50-<?200 differs from PCA pl82-<?200 in the transition 101 -> 
111, that lacks in the former. It means that in PCA pl50-g200, clusters 11- • -1 
cannot coalesce upon encountering. Preliminary numerical simulations indicate 
that this PCA suffers an extinction-survival phase transition by p ~ 0.572, in the 
vicinity of the bump seen in Figure 2. The density of active cells in PCA pl50- 
q200 is smooth at p > p c , without features. A possible scenario for the dynamics 
around p ~ 0.572 in PCA pl82-(7200 is then that the effect of sites surviving 
and spreading together with that of clusters 11- • -1 coalescing increase the density 
of active sites by producing larger clusters that both diminish the erosion from 
the boundaries (by diminishing the total number of boundaries) and create more 
stable segments 111. As p increases from p ~ 0.572, boundaries fluctuate more 
widely, eventually leading to the annihilation of previously existing clusters that 
cannot recover easily, leading to a decrease in the density of active sites and to 
the non-monotonic behaviour seen in Figure 2. Increasing p makes coalescence to 
become more probable and single active cells in triplets 010 more stable, leading 
to a resume of the increase in the density of active sites. 

We shall henceforth focus on the transition around p ~ 0.48 only. 
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4.2 Critical indices 



We determined the critical point p c through an examination of the stationary den- 
sity of active sites pl(p) near p ~ 0.48 for different lattice sizes and also by perusal 
of the scaling relation 

p L (t) ~ r Plu w^{et l l v w, t v± h\/L), (8) 

with e = p — p c ^ [16]. For a very large system, this scaling relation assumes 
the form p(t) ~ t~^ u ^(et 1/v W), with $(x «1)~ constant and $(x > 1) ~ 
or. The first approach gives p c and /3 more directly but is harder computationally, 
due to the very long simulation runtimes to attain stationarity close to the critical 
point. The second approach is more computer-friendly and usually more precise 
but requires first the simultaneous determination of p c and 5 = (5/v\\ and then 
the determination of u». The third exponent u± (or, equivalently, z = v\\fv±) is 
obtained by finite-size scaling. 

Figure 4 shows the estimated critical points p c (L) obtained directly from the 
stationary density profiles pl(p)- A linear regression (LR) fit to the data gives 
p c (oo) = 0.4880(1) (with a puny R 2 = 0.776, though). The inset in the figure 
exhibits the data from the mean-field approximation of orders n = 2, 3, and 4 
(remember that the n = 1 and n = 2 approximations give the same results). The 
extrapolated p c °°^ = 0.5138 from an LR fit to these data (R 2 = 0.984) overshoots 
the more reliable value obtained from the Monte Carlo simulations and is given 
here just to illustrate the slow convergence of the mean-field approximation with n. 

As we have already mentioned, the determination of /3 through a log-log plot 
of p(e) ~ e 13 close to the critical point is a poor choice of method. A better option 
is to plot the time-dependent profile Ph{t) close to p c for some large L. On the 
critical point, p(t) ~ t~ s and we can estimate p c and 5 simultaneously by plotting 
log 6 [pi,(t /b) I ph{t)] against 1/t for some small b. Our data for pi,(t) appear in 
Figure 5a and represent averages over 1000 realizations of the dynamics, all with 
initial density pl(0) = 0.2 in a lattice of L = 16000 sites. From these data 
we estimated p c = 0.48810(5). The data for the instantaneous values of 5 using 
b = 10 appear in Figure 5b. Our best estimate for this exponent comes from the 
curve with p = 0.4881 and is given by 5 = 0.17(1). 

The exponent m can be obtained by plotting t s pi(t) versus te v W and tuning 
to achieve data collapse with different e. Similarly, by plotting t s pi{t) versus 
t/L z atp c for different L and tuning z until data collapse we can obtain z = u»/u±. 
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Figure 4: Critical points p c (L) for several lattice lenghts 8000 < L < 48000. A linear fit 
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Figure 5: (a) Logarithmic plot of p L (t) at fixed e with L = 16000 averaged over 1000 
samples. From the lowermost curve upwards, p = 0.4879, 0.4880, 0.4881, 0.4882, 0.4883, 
and 0.4884. From this dataset we estimated p c = 0.48810(5). (b) Instantaneous values of 6 
obtained from the data depicted in Figure 5a. In panel (b) the uppermost curve corresponds 
to p = 0.4879. From the curve with p = 0.4881 we estimated S = 0.17(1). The dashed 
line corresponds to the best value available for <5dp- 
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Figure 6: (a) Data collapse of the scaled time-dependent density profiles for e = ±0.0001, 
±0.0002, and 0.0003. The upper (lower) branches correspond to p > p c (p < p c ). The best 
data collapse was obtained with p c = 0.48812, 5 — 0.16, and v\\ = 1.70. (b) Finite-size 
data collapse of the scaled time-dependent density profiles on the critical point p c = 0.4881 
for 2000 L < 16000. Best data collapse was obtained with 8 = 0.165 and z = 1.55. 



The collapsed curves for our data appear in Figure 6. Figure 6a was obtained with 
p = 0.48812, 5 = 0.16, and v\\ = 1.7. We can comfortably estimate v\\ only within 
±0.05, so that our best estimate for this exponent is v» = 1.70(5). Combining 5 
and v\\ furnishes /3 = 5u» = 0.29(2). Data collapse in Figure 6b was obtained 
for data taken at p = 0.4881 with 5 = 0.165 and z = 1.55. As before, we could 
determine z only within ±0.05, so that our best estimate for it reads z = 1.55(5). 
The exponents 5, un, and z together determine the universality class of critical 
behaviour of the model, other exponents following from well known hyperscaling 
relations [16]. 

The best values available for the critical exponents of the DP process on the 
(1 + l)D square lattice are 5 DP = 0.159 464(6), ^ DP = 1.733 847(6), /3 D p = 
0.276 486(8), and z DP = 1.580 745(10) [17]. Thus, within the error bars our data 
put the critical behaviour of PCA pl82-g200 in the DP universality class. Our 
estimation of 5 is somewhat high and limitrophe to the accepted 5dp- We probably 
underestimated the uncertainty in this quantity. This higher value of 5 impacted 
the estimation of uu and z, which turned out to be smaller than the accepted values, 
but in these cases still within the error bars. Anyway, we can safely establish the 
universality class of critical behaviour of PCA pl82-<7200 as that of the DP process, 
as no other universality class has exponents close to the values found for it. 
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5 Summary and perspectives 



We showed that the mixture of CA rules 182 and 200 into a PCA suffers an 
extinction-survival-type phase transition at p c = 0.48810(5) in the (1 + 1)D di- 
rected percolation universality class of critical behaviour. Although PCA pl82- 
q200 has two absorbing configurations, 00- • -0 and 11- • -1, this second absorbing 
configuration is hardly achieved, since neither CA 182 nor CA 200 have a station- 
ary density of active sites close to 1 except if the initial configuration is 11- • -1 
itself. We also found that PCA pl82-g200 has some unusual features, notably its 
slow dynamics for very small p, when it approaches the stationary state diffusively 
with a large relaxation time proportional to 1/p, and the bump in its density profile. 

We would like to delve further into the small p behaviour of PCA pl82-g200. 
In this regard, it would not be entirely without interest to consider a stochastic 
reaction-diffusion version of the PCA consisting of the reactions 100 f± 110, 
001 f± 011, and 010 — > 000, making up a sort of constrained, partially reversible 
contact process. In this model, clusters 11- • -1 never coalesce and can only be 
eroded from the boundaries. It is also related with an RSOS interface model in- 
troduced some time ago [18]. We guess that, depending on the rate of the process 
010 — > 000, it displays an extinction-survival phase transition. Moreover, simple 
diffusion 10 f± 01 may represent a relevant perturbation for this lattice gas, since 
it destabilizes 11- • -1 clusters and provides room for more 010 — > 000 reactions. 

There is plenty of room for exploratory incursions into possibly interesting 
single-parameter composite PCA. Of particular interest would be to find composite 
PCA displaying phase transitions in the even branching and annihilating random 
walk or the directed Ising universality classes of critical behaviour [19]. 
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